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ABSTRACT 

We use the IRS on Spitzer to observe the southern part of the reflection nebula NGC 2023, mcluding 
the Southern Ridge, which is a photodissociation region par excellence excited by HD 37903. Five 
pure-rotational H2 emission lines are detected and mapped over and around the Southern Ridge in 
order to compare with predicted level column densities from theoretical PDR models. We find very 
good agreement between PDR model predictions and emission line intensities and ratios measured 
with Spitzer, leading us to conclude that grain photoelectric heating sufficiently warms the gas to 
produce the observed H2 line emission via collisional excitation. On the Southern Ridge, we infer a 
hydrogen nucleus density nn « 2 x 10^ cm~^ and radiation field strength x ~ 10"' relative to the local 
Galactic interstellar radiation field. This high value for x independently predicts a distance toward 
HD 37903 of 300 pc, and is consistent with the most recent HIPPARCOS results. Over the map we 
find that both rtu and x vary by a factor of '^3. Such 2-D variations provide clues about the underlying 
3-D structure of the Southern Ridge field, which appears to be the tip of a molecular cloud. We also 
map variations in excitation temperature and the ortho-to-para ratio, the latter attaining values of 
~1.5 — 2.0 on the Southern Ridge, and find that PDR modeling can readily reproduce observed ortho- 
to-para ratios that are <3 for rotational excitation dominated by collisional processes. Last, the stars 
Sellgren C and G are discovered to be resolved on archival HST images into two point sources each, 
with separations of ^0'.'5. 

Subject headings: infrared: ISM — ISM: clouds — ISM: individual objects (NGC 2023) — ISM: 
molecules — photon-dominated region (PDR) — Stars: individual (HD 37903) 



1. INTRODUCTION 

Photodissociation regions (PDRs) are regions in inter- 
stellar clouds in which far-ultraviolet {6 < hv < 13.6 eV; 
FUV) radiation plays a si gnificant role in th e heating 
and/or chemistry (Ticlcns fc HollenbachI [19851 ) . For ex- 
ample, PDRs are found in reflection nebulae and molec- 
ular cloud surfaces, where the radiation from nearby 
OB stars illuminates the clouds. The incident starlight 
is absorbed by dust and polycyclic aromatic hydrocar- 
bons (PAHs) and is mostly used to excite the PAHs 
and heat the grains. However, a fraction (~0.1 — 1%) 
of the absorbed FUV starlight may be converted into 
energetic photoelectrons, which are ejected from PAHs 
and grains, and subsequently heat the gas. The strong 
FUV radiation acts as a beacon to illuminate the cloud 
structure, and to photodissociate, ionize, and excite gas 
phase chemical species, which otherwise would not be 
seen in emission. Thus, PDRs emit strong far-infrared 
continuum emission from grains, as well as infrared, 
submillimeter, and millimeter-wave line emission aris- 
ing from the warm gas. The FUV radiation can af- 
fect the chemistry in molecular clouds to a depth of 
Ay ^ 5 by maintaining the oxygen that is no t tied 
up in CO in atomic fo rm (Tielcns & HoUcnbach 1 19851 : 
iHollenbach eFall [20091 ). This depth is comparable to 
the mean column density in giant molecular clouds 
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([Solomon et al.l[T987l : [Hever et al.l[2009[ ). The same PDR 
physics that is at work at the surfaces of mo lecular clouds 
also a cts in the diffuse interstellar medium (Wolfire et alj 
I2003f ) and, therefore, much of the ISM is found in PDRs. 



In general, the theo ret ical models 



Sternberg fc Dakarnol [19951: IWolfire et all [20031 : 
Kaufman et al.ll2006HLe Petit et al.ll2006n do a good job 



in predicting the atomic fine-structure line intensities 
and line ratios in PDRs; however, several observations 
of lin e emission frorn H2 usi n g ISO £r immermann et aT] 
[l99l iFuente et al] [l999l: IPraine fc Bertoldi. ,20001) 
and Spitzer ([Goldsmith et al.l 120 lOt ) seem to indicate 
temperatures higher than those predicted by grain 
photoelectric heating alone. In additi on, observ ations 
and modeling of high- J lines of CO (J affe et all [l990l : 
iSteiman-Cameron et al.l [T997D would also indicate that 
models underestirn ate the gas temperature (Tgas). 
iHabart et al.l ([2011[ ) found order-of-magnitude discrep- 
ancies between their PDR model results and Spitzer 
H2 data for rotational level s J > 3 to ward mainly low- 
excitation PDRs. Dedes et al.l ~(|2010D were able to fit 
observations of high- J CO lines by using spherical PDR 
models for an ensemble of clumps dis tributed in size 
and mass. iWeingartner fc Draind ([1999f ) have suggested 
that radiation forces on grains increase the dust/gas 
ratio in PDRs thus leading to enhanced heating rates. 
Dissipa tion of turbulence ([Falgarone et al.l [2005) and 
shocks ([Habart et al.|[20ll[ ) might also be an important 
source of heating in low-FUV field environments such as 
the Taurus molecular cloud ([Goldsmith et al.l l2010t ) or 
in the diffuse ISM. 
Molecular hydrogen, H2 , is a sensitive probe of PDRs 
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in our Galaxy (lAllers et al.ll2005[) a s well as in highly red- 
shifted galaxies ( Sheffer et al.l2009l) . The pure- rotational 
(y = 0) transitions of H2 in the mid-IR are readily ob- 
served by space ins truments, such as the IRS on Spitzer 
(jHouck et all 120041 The high spatial resolution (<10") 
of our Spitzer observations more clearly isolate the emis- 
sion structures and different physical regions compared 
to the previous ISO observations. In general, transi- 
tions between low- J levels of the v — state probe Tgas 
owing to their low critical densities, while the higher- 
excitation energy levels are pumped by FUV radiation 
and probe the radiation field strength and gas density 
(jSternberg fc Dalgarnolll989HBurton et al.lll990D . 

This project has two main goals: first, to use H2 obser- 
vations from Spitzer-IRS of the reflection nebula NGC 
2023 to derive average gas physical conditions in the 
warm molecular regions and, second, to assess the re- 
liability of PDR models to model the gas emission. In 
§ 2 we track the changing values for the estimated dis- 
tance of NGC 2023, present recent imaging results of the 
nebula, and describe our spectral mapping of the South- 
ern Ridge (a.k.a. filament or bar, hereafter SR) area in 
NGC 2023 based on Spitzer data collected by the IRS. 
The intensities of pure-rotational emission lines from H2 
are measured and converted into level column densities. 
These observables are co mpared in 3 with mo del pre- 
dictions based on updated [Kaufman et al.l ()2006[ ) models 
to show that a high level of agreement exists between the- 
ory and observation. Model parametrization allows us to 
map the 2-D distribution of total gas density and FUV 
strength, and to glimpse the 3-D structure of the nebula 
based on model fits. Section 3 also presents maps of the 
ortho- to para-H2 (o- to P-H2) ratio (OPR) and of ex- 
citation temperatures (Tox) for H2, and concludes with 
a discussion of detected fine-structure line emission from 
Si II and Fe II. Section 4 is dedicated to the role of pho- 
toelectric heating compared to other processes. Finally, 
a concluding section provides a textual closure for the 
paper by emphasizing our main results. 



2.1. 



2. THE TARGET AND THE DATA 
What is the Distance to NGC 2023? 



As a reflection nebula surrounding the hot B1.5 V star 
HD 37903, the issue of the distance of NGC 2023 is in- 
timately tied to the distance of its central star, as well 
as being an essential ingredient in estimating the phys- 
ical separation between HD 37903 and the SR. Most of 
the studies dealing with NGC 2023 over the last 30 years 
have employed a consistent range of values of 450 — 500 
pc for the distance toward this object, as well as toward 
the entire nebular complex in Orion. However, assuming 
the depth is roughly the same as the projected extent, the 
angular size of the constellation of Orion predicts a depth 
to distance ratio of '--^35%, amounting to a range o f ^140 
pc for a central value of 400 pc. Indeed, ,Anthonv- Twarogl 
(|1982| ) derived a statistical distance toward Orion B stars 
of 380 ^80° pc. This sho rter distance scale was employed 
by the iWvrowski et al.l (|1997[ ) and iMartini et aD (| 19991 ) 
studies of NGC 2023, but it was quickly abandoned in 
favor of the longer one b ased on the HIPPARC OS par- 
allax of 2.12 ± 1.23 mas (jPerrvman et al.lll997[ ), placing 

HD 37903 at a rather imprecise distance of 470^^70 P*^- 
However, very accurate parallax measurements with 



very long baseline interferometry show the distance to 
the Orion Nebula to be 400 pc, with uncertaint y of 2 — 
6% (jSandstrom et al.l 120071: iMenten et aIll2007D. Addi- 
tional support for a shorter scale is provided bv lCaballerd 
(|2008| ). who derived a distance of 334^22 (or, less likely, 
385 ± 15) pc toward a Ori, a star in the angular vicinity 
of HD 37903. Moreover, the cloud complex that includes 
NGC 2023 and the Horsehead Nebula appears to be in 
front of a Ori (Mookcrjca ct al. 2009). We shall adopt 
350 ± 50 pc as the probable di stance toward HD 37903 as 
suggested bv lMookeriea et al.l (p0'09i) . Clearly, HD 37903 
is located at the very near side of the Orion nebular com- 
plejfl, implying that the SR is separated from the star 
by a projected distance of (4.0 ± 0.6) xlO^"^ cm, or 0.13 
± 0.02 pc. 

2.2. Non-Spitzer Images 

In order to become familiar with the appearance and 
structure of NGC 2023, we include here two high-quality 
public data products that have been obtained in recent 
years. The first is a near-IR image obtained at ESO (left 
panel of Figure [T]) , which clearly shows the structure of 
the nebula that heretofore had been hidden from visual 
view. The star HD 37903 has carved a quasi-spherical 
cavity into the dense molecular material out of which 
it formed. Consequently, its energetic FUV radiation 
illuminates the ridges of high-density gas, producing a 
reflection nebula and emission ridges that are detected in 
H2 line and dust (as well as PAH) continuum emission. 
Sellgren D, the eye-catching orange object below the pink 
SR, is resolved into at least three sources, whereas it was 
seen as an elongated object in previous renditions. 

A portion of an archival image from the Hubble Space 
Telescope is also shown in Figure [TJ Taken in red visible 
light with the ACS, it is the highest-resolution image ever 
taken of NGC 2023. The SR is resolved into an ensemble 
of clumps, all part of massive ridges of cloud tops that are 
basking in the FUV starshine of HD 37903. The SR was 
already known fro m ground-based work to be ^2" wide 
()Field et al.lll994 (1998). but the newly resolved clumps 
are smaller than that by a factor of a few. They are, on 
the other hand, much larger than the O'.'OS angular size 
of each pixel in the ACS field of view. 

As a bonus, the (full) HST image also resolves two 
stars into their constituent components for the first time. 
These are Sellgren C near the western end of the SR, and 
Sellgren G, which is labeled on the ESO image. How- 
ever, we cannot be certain that two actual binaries are 
involved, since young stellar objects, such as these PMS 
stars (Sellgren 1983), may be surrounded by a residual 
edge-on debris disks from the formation p rocess, which 
might mimic a binary source. According to lDePoy et al.l 
()1990( ). their star 9 (Sellgren G) is consistent at the l-cr 
level with standard reddening vectors, i.e., it is not a con- 
vincing case of special extinction by a PMS "cocoon" . In 
both cases the components are 0'.'4 ± O'.'l apart, or pro- 
jected separations of ^150 AU. 

2.3. Spitzer Spectroscopy 

We acquired Spitzer-IRS spectra of H2 emission lines 
toward NGC 2023 by employing the three modules SL, 

The short distance scale has now been confirmed by HIPPAR- 
COS, see end of 5 5 for details. 
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SH, and LH, as listed in Table 1, where the second letter 
'L' or 'H' stands for low- or high-spectral resolution, re- 
spectively. The lower-resolution SL module is split into 
orders SLl and SL2 and provides potential coverage of 
H2 emission lines from J = 4 — 9 levels, i.e., of transitions 
S(2) through S(7). However, owing to substantial blend- 
ing with strong PAH features, we could detect only S(5) 
at 6.91 /zm in SL2 data, and S(3) and S(2) at 9.66 and 
12.28 /zm, respectively, in SLl spectra. SH also covers 
S(2), but at resolution wlO times higher than that of SL, 
and has exclusive coverage of S(l) from J = 3 at 17.03 
/zm, whereas LH provides the only coverage of S(0) from 
J = 2 at 28.22 yitm. 

The data were processed into cubes using the match- 
ing Spitze r Science Center pipeline version (S18.7.0) and 
CUBISM (jSmith et al. 2007) version 1.7. AU non-LH ex- 
posures were spatially degraded by re-gridding onto the 
(larger) pixels of the LH field of view prior to further 
analysis. Figure [2] shows the proper celestial location 
and orientation of the LH field over an 8 /xm image from 
Spitzer-IRAC. With 15 x 14 pixels in the LH field, each 
4'.'46 across, the irregular border of the combined inter- 
section area of all modules measures w 1' x 1' along the 
instrumental X and Y directions. 

Emission line maps were constructed by using 
IDL/GAUSSFIT to fit line profiles and to derive inte- 
grated line intensities that included continuum fitting 
and removal. Table 2 lists derived GAUSSFIT param- 
eters and their uncertainties for the case of map medians 
(dominated by off-SR pixels) and for the case of the single 
on-SR pixel LH[7:8], which coincides with the location of 
peak H2 emission. Figure [3] shows the spatial distribu- 
tion of H2 emission for the five detected transitions S(0), 
S(l), S(2), S(3), and S(5). In these maps, the SR is a very 
prominent source of H2 line emission, with a partial cov- 
erage of additional emission from the South-Southeastern 
Ridge toward the instrumental top right corner. (Note 
that here, and in following Figures, the instrumental ori- 
entation of the LH map is to be employed in order to 
avoid both additional interpolation of the data and the 
wasteful white margins inherent in celestial orientation.) 
The H2 shows good agreement in position and orienta- 
tion with the SR seen in the IRAC 8/im image, which 
traces mainly PAH emission (Figure [2]) . 

The target area observed by Spitzer toward NGC 2023 
shows prominent and broad emission features from PAH 
molecules together with strong and narrow (unresolved) 
emission lines from H2, see Figure IH This spectrum 
was obtained by averaging 15 LH pixels that sample 
the emission from the SR in all four IRS modules. Two 
minor contributions from atomic species on the SR be- 
long to [Fe II] and [Si II], at 25.99 and 34.82 fim, re- 
spectively, which are expected to be PDR observables 
(jTielens fc HoUenbachl 119851 : Kaufman etlUI [2006,) and 
are presented in § 3.4. Two other fine-structure transi- 
tions commonly found in H II spectra, [Ne II] at 12.81 
fim and [S III] at 18.71 /im, are detected at extremely 
weak levels of emission all over the map. We shall con- 
centrate in this analysis on the H2 lines, comparing them 
with model predictions. The analysis of the PAH features 
shall be presented in a later paper (Peeters, et al. 2011, 
in preparation). 

2.4. IRS Calibration Issues 



2.4.1. Differential Extinction 

We applied extinction corrections to our data consis- 
tent w ith previous st udies of H2 line emission from NGC 
2023. iBurtonI (jl993i ) found that Ak = 0.3 mag from 
their vibrationally excited emission lines of H2, yield- 
ing Ay = 2.8 — 2.3 ma g for the range Ry = 3.1 — 5.5 
()Mathisl[l990l ). whereas iBurton et all (|1998D mentioned 
that Ay is likely to be ~3 — 5 m ag, or Ak ~ 0.3 — 0.65 
mag. iDraine fc Bertoldil (|1996[) and iDraine fc Bertoldil 
(|2000( ) adopted Ak of 0.2 and 0.5 mag, respectively. For 
our corrections we adopted Ak = 0.5 mag, which corre- 
sponds to Ay = 3.8 — 4.7 mag; 

Using the extinction curve of lMathii (|199O0 , the extinc- 
tion correction varies between 6 — 29%, with the largest 
correction applying to S(3). The latter value is simi- 
lar in magnitude to systematic effects in the absolute 
calib ration of Spitzer-IRS fluxes, estimated to be ~20 — 
25% (jGalhano et al.l[2008l : IDale et al.l[2009h . For the ro- 
tational levels J = 2 and 5 with the lowest and high- 
est extinction corrections, respectively, we determined 
that varying both Nj values by ±20% shifts their de- 
rived T52 value (i.e., Tex of level J — 5 relative to level J 
= 2) by -1-3 and —8%. This smaller change, owing to a 
weak dependence on the ratio of colu mn densitie s (sinc e 
T52 oc ln[iV2/A^5]), is consistent with IDale et all (|2009| ). 
who reported ±10% for derived line ratios extracted from 
Spitzer-IRS data. 

2.4.2. Background Subtraction 

Background contribution from zodiacal emission is al- 
ways present toward celestial targets, as a function of 
space (direction) and time (date). Toward NGC 2023, 
dedicated background exposures were taken only for the 
SL observations, revealing a wavelength-averaged ratio 
of background to data of 10%. We employed the zodi- 
acal emission calculator from SPOT, for each date and 
pointing center of each observation in order to compute 
the expected level of background emission for our data, 
and found good agreement between observed and calcu- 
lated zodiacal emission for the SL data. Thus we are 
assured that subtraction of calculated zodiacal emission 
from data that lack background exposures does not lead 
to any significant errors. Furthermore, since our analy- 
sis concerns emission intensity following continuum sub- 
traction, such intensities are insensitive to the presence 
of continuum-like background levels. 

2.4.3. Intensity Inter- Calibration 

There is a significant wavelength overlap at 9.97 — 
14.74 /im between the SH and SLl modules. Previous 
works have shown a continuum mismatch between mod- 
ules, with different modules requiring different scale fac- 
tors to_briiig_tlie continuum into agreement. For exam- 
ple, iBrandl et all (|2004). had to shift their SH data by 
+36% and their SLl data by +17% relative to LL data, 
showing that in their case, the SLl scale was higher by 
+16% than the scale of the SH module. 'Quan z et al.l 
(2007), on the other hand, found higher readings from 
SH data relative to SLl data, and attributed these 8 — 
25% scale shifts to different sht orien tations relative to 
source emission. iBeirao et all pOOSi ) found up to 50% 
differences in fits of PAH features based on SH and SLl 
spectra. 
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We performed a pixel-wise analysis of SH data ver- 
sus SLl data, which included background subtraction as 
well as SH resolution degradation to the lower spectral 
resolution of SLl data. The absolute intensity scale of 
the SH data was found to be higher by 19% than that 
of SLl data. Since we consider the higher signal-to-noise 
SLl data to be more reliable spectrophotometrically, the 
SH (and LH) data were re-scaled by 0.84 prior to mea- 
suring integrated intensity values for the S(0), S(l), and 
S(2) lines. 

The combined uncertainty in our measurements owing 
to the dominant calibration issues of absolute flux uncer- 
tainty and inter- module uncertainty is <30%. 

2.5. Comparison of Spitzer IRS with ISO SWS 

The extinction-corrected H2 line intensities were 
converted into level column densities using Nj = 
AttIj/AjAEj cm~^, where Nj,Ij,Aj, and AEj stand 
for the column density, emission intensity, Einstein A- 
coefRcient, and transition energy for each upper level J. 
The H2 rotational emission lines are quadrupole transi- 
tions and thus an optically thin conversion is appropriate 
for all environments in the ISM. For mapping purposes, 
this procedure was followed for each pixel over the field 
of view. To compare with previous ISO observations, we 
employ the mean on-SR H2 column density based on Ij 
values from 15 pixels that sample SR emission. These 
pixels were also employed in the extraction of the aver- 
age on-SR spectrum shown in FigurelH The SR-averaged 
log Nj values for J = 2, 3, 4, 5, and 7 are 20.29, 19.77, 
18.95, 18.60, and 17.48 cm'^, respectively. 

/5'0-SWS observations of pure rotational transitions 
toward the SR of NGC 2 023 were shown in Figure 5 of 
iDraine fc Bertoldil (|2000l ). from which we extracted Nj 
values for the five emission lines detected in our Spitzer 
data. These values were then rescaled by 1/1.8 in or- 
der toj;emove_the arbitrary beam filling factor employed 
bv IDraine fc Bertoldi (2000), resulting in log Nj values 
of 20.44, 19.84, 18.90, 18.56, and 17.24 cm'^. Conse- 
quently we find ISO to Spitzer column density ratios of 
1.4, 1.2, 0.9, 0.9, and 0.6, averaging 1.0 ± 0.3, which has 
the uncertainty expected for two data sets with ^20% 
uncertainty each. Still, it is fascinating that the ratio 
appears to have a monotonic decline with increasing J, 
possibly related to differences in beam sizes employed for 
d iffering wavele ngt h regi mes on both spacecraft. 

iHabart et al.l ()2004D quote the (extinction- 
uncorrected) ISO observable / — 1.65 x 10~^ erg 
s~^ cm~^ sr~^ for the S(3) line, taken with a 19" 
beam centered on a point 60" due south of HD 37903. 
This position lies north of the SR and prevents us 
from directly comparing it with SR intensity values, 
however it still lies within our mapped region. We 
simulated this beam on our IRS map and found that 
the Spitzer H 2 emission is 4 2% higher than the value 
given by Hab art et al.l (j2004f ). hence a column density 
ratio of ISO / Spitzer = 0.7, which is consistent with the 
ratios given above. The agreement between IRS and 
SWS results is quite good considering the uncertainty 
in the ISO beam filling factors and extended source 
calibrations for both ISO and Spitzer. 

3. MODELING SPITZER DATA WITH PDR MODELS 



3.1. Model Parameters 

In order to test the results of PDR models, we compare 
-^j(H2) values derived from our observed emission lines 
toward NGC 2023 with values p redicted by our PDR 
model ()Kaufman et al.lll999l[2006l ). This model includes 
the calculation of H2 processes from the Le Petit et al.l 
((2OO6) code as discussed in Kauf man et al.l ([20061 ) . These 
include radiative excitation and dissociation, dissociation 
heating, collisional excitation and deexcitation, radiative 
cooling, and heating by deexcitation of excited levels. 
Th e 0-H2 to P-H 2 conv ersion on grains is included in 
thelLe Petit et al.l (|2006f ) code as described in lLe BourlotI 
(12000). We also include a factor of 2 enhanced H2 for- 
mation rate as sug gested bvlHabart et al.l (|2004| ) and dis- 
cussed in lKaufman et al.l (I2006D. The OH and C O chem- 
istry has been updated as in lWolfire et al] (|2010i ) and we 
test models for normally-incident photons. The model 
output consists of the H2 column densities in each J 
level integrated along the normal to the PDR surface, as 
well as the normally emitted H2 line intensities. The ro- 
vibrational quadrupole transitions are all optically thin. 
We will vary two main model parameters to obtain a best 
fit to the observations while holding all others constant. 
These are the hydrogen nucleus density, n-g, and t he in - 
cident radiation field, Xi in units of the iDraind ()1978f ) 
interstellar radiation field. We use the notation that x 
is the ratio of FUV field incident on the surface of the 
PDR divided by the free-space field in the local interstel- 
lar radiation field. Thus x — Fpuy /AttIj:), where Fpuv 
is the incident FUV flux and Id = 2.2 x 10"'* erg s~* 
cm~^ sr~^ is the Draine intensity for the local ISM. 

In order to establish initial parameter values for nn and 
X^ as well limit their variation during modeling to values 
that are consistent with known ranges of other observ- 
ables, we need independent m ethods for estimating nu 
and X- iWvrowski et al.l ()2000l ) applied PDR model re- 
sults to their observed angular separations between emis- 
sion from H2 and C91a toward NGC 2023. They found 
a range of nn = 0.6 — 1.4 xlO^ cm""^, which agrees with 
previous findings that this is a relatively high-density 
PDR and provides us with an initial density value of 
nn « 10^ cm-3. The FUV field reaching the PDR gas 
from HD 37903 depends on the latter's luminosity and 
the projected (lower limit) physical separation between 
the two, which depends on the distance of NGC 2023 
from Earth. For a B1.5 V star of 12 Mq (|Conti et all 
12008") we find an FU V luminosity of 1.13 x lO^L©, based 
on Parrav ano et al.l ()2003 1 . An angular separation of 78" 
between the star and the SR results in an FUV field 
strength of x = 9.27 x 10^/D^, where D is the distance 
in pc to HD 37903. Using the adopted value of D = 350 
± 50 pc (§ 2.1) predicts a range of fairly high values for 
X of 7.6l^'g X 10"^. We shall explore a range of x val- 
ues corresponding to ^2-cr variation around this central 
value. 

In comparing normal model line intensities to obser- 
vations there are several factors to keep in mind, which 
can be divided into two classes. One class includes fac- 
tors that may increase the ratio of observed intensity to 
model intensity. The first factor is the possible presence 
of multiple PDRs along the line of sight, /p. A second 
factor is the inclination angle {9) of a single PDR layer 
relative to the line of sight. Limb brightening is expected 
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to occur for any PDR inclination 6* > 0°, causing the 
source line intensity to appear brighter than a face-on 
[e = 0°) model PDR by the factor fg = l/cos{e). Both 
factors can be lumped together as a single raising factor, 

/+ = hfe > 1. 

The other class to consider includes factors that may 
decrease the ratio of observed intensity to model inten- 
sity. First is the fraction of the beam area that is filled 
by the source emission, /b- A second factor involves the 
incidence angle of FUV illumination, which affects the 
penetration depth of the radiation field. The intrinsic 
(deprojected) geometry of the exciting star (HD 37903) 
to the SR should produce a range of values between 
strictly normal (0 = 0°) and strictly parallel (0 = 90°) 
incidence, depending on the details of the shapes of gas 
clumps on the SR. For any value of > 0°, the in- 
tensity of observed H2 emission is cut down by cos{(f>) 
owing to the oblique path of the FUV photons through 
the gas layer, and thus we have for the lowering factor 

/- = /b/0 < 1. 

The effects of the raising and lowering factors obviously 
work in opposite directions, such that the total effective 
ratio of data to model is fcs = f+f- = /p/e/B/0. Effec- 
tive ratios for map pixels will be derived as logarithmic 
differences (or model "shifts") in § 3.2 via matching of 
absolute model Nj values between models and data. The 
value of fcs only tells us whether /+ or is the dom- 
inant effect, but not the actual values of any individual 
factors involved. 

3.2. Absolute Column Densities as a Function of J 
3.2.1. A Single Normal Model 

The absolute column density values, Nj{ll2), consti- 
tute our primary means of comparing PDR model results 
with the observations. Figure [5] illustrates such a com- 
parison for a single on-SR pixel, where a run of observed 
Nj values vs Ej, the level energy above the ground state, 
is compared with output from the best matching model. 
For illustration purposes, data are shown both before and 
after extinction correction is applied. The best match be- 
tween models and data was determined from the smallest 
root mean square deviation (RMSD) of the differences in 
dex between modeled and observed absolute Nj values. 
In this way, both relative Nj ratios between even- and 
odd-J levels, as well as global absolute Nj values, were 
being fitted with the best-matching model for each pixel. 
The RMSD-minimizing search for fcs was performed via 
globally shifting the Nj values of each (tih, x) model in 
steps of 0.01 dex relative to the data, over a range of ±1.0 
dex. (Larger ranges of up to ±10 dex were tested, but did 
not improve the fits.) With a minimum RMSD value of 
0.050 dex (±12%) in FigureO it is obvious that very good 
agreements can be found between observed and modeled 
column densities in terms of the overall shape of the Nj 
curve. 

As a test case, we fit several regions with fixed values 
of riH = 10^ cm-3 and X = 5 X 10^. The median RMSD 
value on the SR is 0.144 dex, or a difference of 39% be- 
tween data and models, i.e., clearly not fitting the data 
within the expected observational uncertainties. See Ta- 
ble 3 for a listing of results for fixed tih and x that include 
other regions around the SR. Owing to the poor fit by 
a single model of fixed parameters, better agreement be- 



tween data and models is expected from expanded ranges 
of the two parameters riu and x- 

3.2.2. Expansion to a Multi-Model Grid 

In order to reveal density and FUV flux variations, 
and thus take advantage of the mapping performed by 
the IRS, we have generated a grid of models over the nn 
vs X parameter space. The grid of 49 normal models 
used here covers 0.6 dex in nu and 0.6 dex in x with 
a step size of 0.1 dex. It is a subset of a larger grid 
comprising of 210 models, most of which did not fit any 
of the observed Nj curves. The parameter ranges of nn 
= 0.5 — 2 xlO^ cm~^ and x = 4 — 16 xlO'^ were found 
to provide much better fits (smaller RMSD) than using 
a single model over the entire map. The results of this 
multi- model fit are shown in Figure IHl 

In this Figure we see that the highest density values 
of 2 X 10^ cm~'^ are found on the SR, as well as on its 
neighbor, the SSE ridge. The range of nn over the en- 
tire map is in v e ry go od agreement with that given by 
IWvrowski et al.l (|2000t ). The x map shows an enhance- 
ment of the FUV field on the SR of x = 10^, which is 
~30% higher than 7.6 x 10'^, the central value expected 
according to the distance to NGC 2023 (§ 3.1), but is 
nonetheless consistent with the upper l-cr value for the 
expected x- 0^ fact, our initial comparison with ex- 
pected X values based on the longer distance scale of 
450—500 pc for NGC 2023 had resuUed in model-to- 
expectation ratio of 2.4;']q'2. Such a very significant dis- 
crepancy prompted us to investigate the issue of the dis- 
tance to HD 37903 more thoroughly, as reported in § 
2.1). Taken at face value, a prediction of x = lO"' means 
that HD 37903 could be even closer by 1 sigma than the 
adopted short distance scale, i.e., 300, instead of 350, pc 
awajQ. 

Both parameters diminish farther from the SR, except 
for regions with high x values farther to the North of 
the SR (lower right corner) that have smaller (projected) 
distance from the exciting star, HD 37903. The median 
RMSD value for the entire map area is 0.062 (Table 3). 
Note that for the SR this value has improved by a factor 
of >2 relative to the single-model fit. The highest values 
of the global shift are also found on the SR, where the 
median is —0.13 dex, or /eg = 0.74. Most of the mapped 
(off-SR) area has a value of /eg ^ 0.49, hence 1.5 times 
smaller than the on-SR value. 

3.2.3. Comparing fed wtth Previous Results 

The /cff-corrected RMSD median of 0.062 shows that 
the combination of spatially higher-resolution data from 
Spitzer-IRS and recent improvements in PDR models ap- 
pears to yield a very good level of agreement of 15% be- 
tween observations and predictions of H2 emission. Over 
the map (Fig. ^ we find that -0.57 < log f^s < -0.09, 
or 0.27 < /+/- < 0.81, a variation by a factor of 3. 
Since /+/- < 1 for all pixels, /_ is the dominant factor 
for comparing models to observations of this region (i.e., 

/+ < I//-). , , 

In their modeling of t he SR. [ Draine fc B crtoldi (19961) 
corrected the data from IHasegawa et al.l (|1987. ') by using 

^ This distance of 300 pc has now been confirmed by HIP PAR- 
COS, see end of § 5 for details. 
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a beam filling factor of /b = 1/6, and employed an incli- 
nation correction of fe — 5, oi 6 Ki 78°. Their fcs = /e/e 
= 0.83 is closely matching the higher values found on our 
map. [D raine & Bcrtoldi (2000) used the same fg from 
IDraine fc Bertoldil 11996)), but /b = 1/1.8 = 0.56, result- 
ing in /eff = 2.8, a v a.lue >3 times higher th an our map's 
largest value of 0.81. iKaufman et al.l ()2006[) used a proxy 
for the beam factor in the form of two-phased H2 medium 
(jSteiman-Cameron et al.lllQQTf) . in which the dense gas 
from the SR was contributing at a level of 20%, and thus 
in effect using /b — 0.2. In addition, they used /e = 6 in 
order to achieve a match between data and models. In 
short, this amounts to using fcs = /s/b = 1-2, a value 
higher by ~50% than our map's highest value. The con- 
tributions of multiple PDRs along the line of sight, /p, 
and the angle of FUV illumination, f^, were not explic- 
itly mentioned by these studies. 

We first consider the effects of varying /b while hold- 
ing the other factors fixed. Adopting the assumption 
of .fp/tA = 1, and usin g fe — 5 on the SR, as per 
IDraine fc Bertoldil (fT996l) . we find /b = fcs/fe = 0.74/5. 
= 0.148. However, off the SR, where values of fcs are 
smaller, i.e., fcs = fefs — 0.49 (Table 3), the implied 
beam filling factor would also be smaller, /b = 0.49/5. 
= 0.098. This trend is not likely since away from the 
SR, the source is more extended and not as concentrated 
as the SR, and therefore we expect /b to become larger 
(approaching 1), not smaller. Thus we conclude that it 
is not /b that is driving down the values of /- (and thus 
of /off) away from the SR. In the next section we explore 
the possible 3-D shape of the SR and show how a simul- 
taneous reduction in both and in fg may explain the 
smaller off-SR values of fes- 

3.2 A. Clues for 3-D Nebular Structure? 

Our map exhibits fcs below unity, with values decreas- 
ing away from the SR and possibly indicating a dominant 
reduction in any of the three factors /p, fg, or Z^. AU 
three factors are related to the 3-D structure of the neb- 
ula and their variability over the map should be helpful 
in deciphering the relative configuration between on-SR 
and off-SR regions. 

Let us define a as the angle at the source between our 
line of sight and the direction of incident FUV radiation 
from HD 37903. A contour plot of the product fgf^ = 
cos(0)/cos(6') is shown in Figure [71 based on the rela- 
tionship a — 9 ± (p, or (j) — \a — 6\. The lower half of the 
figure includes the values fgf^ < 1, although geometric 
configurations toward the lower right corner are excluded 
owing to the FUV source dipping below the horizon of 
the FDR face. There are two lines on the map where 
(j) = 0, leading to fgf^ = 1. One is where (f) = 9 — a/2., 
as indicated by the diagonal line, and the other is where 
a = 0°, hence along the ordinate. Between these two 
lines and for any given 9, fgfcf, reaches a maximum when 
0° < a < 90°, along the line = 0° where cos((/)) = 1 
and 9 = a. 

The impression from Figure [1] is that HD 37903 is lo- 
cated at the center of a bubble, and that the SR (as well 
as other ridges of intense emission) is marking the irradi- 
ated boundaries of the bubble. We may thus assume that 
a is not too far from 90°, and thus the value of fgf^ is 
expected to be near its maximal value. Furthermore, fg 
probably achieves its highest value on the SR, assuming a 



near edge-on view of the bubble's boundary. Away from 
emission ridges no such narrow features are seen, and 
thus moving away from the SR would result in lowering of 
the value of fgfcj, — > 1, as (j) 9. For the spherical bubble 
scenario, we expect that 45° ^ a ^ 135°, which results 
in 0.5 < /0 < 1.0 and thus yielding fgf^ > 1 for any 
fg > 2. This constraint is cert ainly true for 9 = 78°, th e 
nominally adopted value from IDraine fc Bertoldil ()1996[) . 

One way to achieve a 3-D configuration such that feS 
becomes smaller with distance from the SR is to assume 
that the molecular cloud banks have a quasi-pyramidal 
cross section. A cartoon that presents such a configu- 
ration is shown in Figure [51 [An earlier, albeit more 
simpli fied, cartoon was given in Figure 4 of [Field et al.l 
(|1994f) .] Thus the southern slope of the SR may face our 
line of sight (with smaller 9 and fg — >■ 1), while the sur- 
face normal is at a larger (/^ — > 0) from the direction 
to HD 37903. In this case, starshine can still hit directly 
in a normal fashion the top of the SR that rise toward 
the star (as can be seen in Figure [T]), but the off-SR gas 
is facing away from HD 37903 and is receiving slanted, 
and thus reduced, levels of FUV flux. In particular, still 
keeping a range of 45° ^ a < 135°, but assuming a 45° 
turn of the viewed slab toward our line of sight (e.g., 
61 = 78° - 45° = 33°) allows for a range of > 12° such 
that 0.00 < f(j, < 0.98. In other words, relative to on-SR 
values, the off-SR pyramidal slope provides a combina- 
tion of smaller 9 and larger that can readily achieve 
fefij> ^ Ij thus accounting for the observed reduction in 
/cff away from the SR. 

In this picture the SR is a manifestation of one ridge of 
clouds forming the outer layer of the much larger, denser, 
and darker molecular cloud. This cloud was outlined by 
HCO+ via the millimeter observations of lWvrowski et aP 
dlOOO), whose Fi gures 3 and 4 clearly show that the FDR 
(H2 and C recombination line emission) is located on the 
HD 37903-facing side of the denser structure. Indeed, 
the more opaque off-SR sight line may also explain any 
potential reduction in /p relative to the on-SR sight line, 
thus providing a 3-D scenario that simultaneously re- 
duces all three factors affecting /off. 

3.3. Excitation Diagrams of H2 
3.3.1. Construction and Interpretation 

While it h as been un derstood since the earliest ob- 
servations (G atlev fc Ka ifu 1987; Hascgawa ct al. 19871) 
that ro-vibrational emission in NGC 2023 is dominated 
by FUV fluorescence, the situation regarding rotational 
emission from the v = Q state could depend on both col- 
lisional and radiative excitations. From our models we 
flnd that the S(0), S(l), and S(2) lines are produced at 
FDR depths where H2-H2 collisions thermalize level pop- 
ulations owing to collisional deexcitation rates that are 
much higher than radiative decays. For the lower Tgas 
involved (<300 K), relevant critical densities of the J — 
2, 3, and 4 levels for H2-H2 collisions are 20, 360, and 
5100 cm~'^, respectively, clearly below the inferred gas 
density of ^ 10^ cm~'^. The S(3) and S(5) lines arise at 
shallower depths near the PDR surface, where the gas 
is warmer and H2-H collisions dominate the excitation. 
With respective critical densities for collisions with H (at 
~600 K) of 3 X 10"* and 1.8 x 10^ cm'^, the upper levels 
J = 5 and 7 are fully and marginally thermalized, re- 
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spectively, by such collisions. This establishes that the 
excitation temperatures of the H2 derived directly from 
observations are more closely associated with the gas ki- 
netic temperature achieved through photoelectric heat- 
ing than with the FUV pumping rate of the H2. 

In addition to (model-derived) gas density and FUV 
intensity, two oth er facets of Hg micr ophysics are the o- 
H2 to P-H2 ratio ([Burton et al.lll99^ and the excitation 
temperatures of rotational-level populations. Values of 
the OPR, Tex, and ground-state population (No) may be 
directly extracted from excitation diagrams employing 
In Nj/gj vs Ej, as shown in Figure HI For gas at LTE 
and with a single value of Tqx, the populations obey an 
exponential (Boltzmann) distribution so that the excita- 
tion function is a straight line having a slope of — I/Tqx. 
However, when the emitting region includes a range of 
H2 excitation temperatures the resulting function will be 
curved. 

Each excitation curve may be approximated by a sum 
of two asymptotic straight lines representing two inde- 
pendent single- Tex H2 components. One is a "cold" com- 
ponent with T™'^ w T32 and the other is "hot" with 
jihot ^ Extrapolation through Ei and Eq provides 
estimates of the ground state column densities A^i and 
Nq. However, this method of estimating the H2 column 
is sensitive to only those regions that are warm enough 
to emit in the rotational transitions and thus neglects 
the colder interior gas. In addition, owing to the mono- 
tonically declining character of the — l/Tex curve, the 
straight line fits tend to overestimate Tgas (via Tip) in 
the emitting regions. 

In the case of H2, the statistical weights gj include a 
2 J+1 factor from orbital statistics and the nuclear factor 
of 3 for odd- J levels. An excitation curve with a misalign- 
ment (zigzag) between odd- and even- J populations is 
the indication that the OPR 7^ 3. The average OPR be- 
tween observed columns in odd- and even- J states may 
be determined by de-zigzagging the curve. A fitting pro- 
cedure was accomplished by shifting the odd- J N values 
toward the even- J levels, with the resulting de-zigzagged 
curve of smallest RMSD providing the observed value of 
the OPR. 

3.3.2. Mapping the OPR 

Under local thermal equilibrium (LTE) and T > 
300 K, 0-H2 is three times more abundant than p- 
H2, hence OPR = 3, owing to nuclear statistics. The 
OPR, which is initially set at the time H2 is formed 
on a surface of a dust grain, can be changed over time 
via collisions between H2 molecules and other gas con- 
stituents including H and H+, and through accretion 
and subsequent ejection from grai ns (Burton ct al. 1992; 
ISternberg fc Neufeldl [1991 iLe Bourloti 12000) . Radiative 
processes in H2 cannot change an original OPR value 
since they involve quadrupole transitions of A J = 0, ±2 
in the ground electronic state. The first evidence for H2 
molecules with OPR 7 ^ 3 in the ISM was provided by 
iHasegawa et al.l ([19871 ). who observed vibrationally ex- 
cited H2 with OPR = 1.4— 2 tow ard the SR of NGC 
2023. ISternberg fc Neufeldl ([19991 ) calculated that an 
OPR value of 1.7 would be observed in vibrationally ex- 
cited H2 even from a gas with 0-H2/P-H2 — 3, owing to 
differences in pumping rates caused by the propensity of 
0-H2 to self shield before P-H2. Thus, definitive evidence 



of OPR other than 3 in a PDR has remained elusive. 
Recent measurements of the OPR in shocked gas, using 
pure rotational transitions of H2 , show more definitively 
that the ratio can differ from 3 owin g to time dependent 
effects of the OPR conversion (e.g., Neufeld et al.l [20061 : 
Mar et et al.l 120091 : [Neufeld et al. 200!^^ 

A de-zigzagging procedure was repeated for all data 
pixels in the field of view, yielding the OPR map in Fig- 
ure [T0| First, the global range of OPR variation over 
the map is 0.9 — 2.0, or wl.S — 2.0 inside the 75% inten- 
sity contour on the SR. This dep arture from OPR = 3 
is very similar to previous results (I Hasegawa et al.l 119871 : 
lBurton|[l99l iFleming et al.ll2010D . Second, the highest 
OPR values are not aligned with SR intensity contours, 
but are found to be shifted toward the north, in the di- 
rection of the exciting star, HD 37903. This phenomenon 
is corroborated by OPR values on the star-facing side of 
the neighboring South-Southeastern Ridge. 

Modeled OPR values are determined by a variety of 
conversion processes between o- and P-H2, with the dom- 
inant ones shown in Figure [TlJ Local OPR values are >3 
throughout shallow PDR layers with Ay ^ 1, but rapidly 
decline to ^1 for Ay > 2. We performed an identi- 
cal de-zigzagging analysis of Nj values from the PDR 
model output. The integrated OPR value was found to 
be 1.8, in agreement with observations. In other words, 
the PDR model, with its various routes of 0-H2 ^ P-H2 , 
provides a close match to the run of Tqx and H2 level 
abundances that closely duplicates the observed OPR in 
the J = 2 — 7 states. This result is important in showing 
that a steady-state PDR model can naturally (i.e., with- 
out adjustments) reproduce observed OPR values that 
are ^ 3, based on integrated columns through the entire 
PDR, without recourse to explanations involving time 
dependent effects. We emphasize that owing to colli- 
sional domination, the H2 rotational populations studied 
here sh ould not be affected by dif ferential self shielding 
effects (ISternberg fc Neufeldl [l999l) . 

3.3.3. Mapping Tex and No 

Values of Tqx derived from our Spitzer data are shown 
in the upper panels of Figure [l2j Comparing the OPR 
map of Figure [TOl with the T^^^'^ map, we see that the 
highest T™''^ values are found on the side of the SR facing 
away from HD 37903, where OPR values are lower, and 
vice versa. 

Overall, the observed range of OPR values is lower 
than their LTE values for the observed range of T™'*^. 
Specifically, a range of OPR = 0.9 — 2.0 may be obtained 
under LTE from a temperature range of 74 — 118 K. All 
observed Tex values derived from these data, as well as 
the Tgas model values (§ 4), are >144 K and thus would 
predict an LTE OPR of >2.4 over the entire map. In 
other words, the H2 has a lower OPR than values cor- 
responding to all indicators of Tgas. (However, as previ- 
ously remarked, extrapolations of T32 to lower- J levels 
tend to overestimate Tio, which is not observed via ro- 
tational lines.) Although this result is derived strictly 
from observed Nj values, we may use the model to ex- 
plain this effect thanks to the very good reproduction 
of Nj, Tex, and the OPR by our PDR model. We sug- 
gest that when Tex and the OPR are calculated from the 
integrated columns, Tex is weighted toward Tgas in the 
line emitting regions, while the OPR is being decoupled 



8 



Sheffer et al. 



from Tgas owing to increasing influence of T^ust- Indeed, 
beyond Ay w 1.5 the OPR is dominated by the 0-H2 — 
P-H2 conversion process via H2 accretion onto cold dust 
grains (Figure [TT]) . These grains have Tdust <C Tgas, thus 
inducing lower integrated OPR values than LTE predic- 
tions based on Tgas- 

The two lower panels of Figure [T^] present the derived 
No populations of "cold" and "hot" H2, as defined by 
the two upper Tex panels. It is seen that colder H2 is 
uniformly spread over the map, but that the hotter H2 
population is found concentrated on the bright ridges of 
highest intensity. This picture of correlation between re- 
gions of higher x values and hotter H2 is confirmed by the 
distribution of T^°* in the same Figure. Thus observables 
and derived characteristics that show the distribution of 
hotter, excited II2 levels, also track regions with higher 
FUV radiation fields. 

From a given value of iVo and a single- Tox level popu- 
lation distribution, the total H2 column density can be 
estimated. Using the partition function formula from 
IHerbst et al.l (Il996l ) gives iV™''! « 0.02^7 N;j°^'^T^°^'^ = 
2.4 X 10^"'^ cm~^. We may then employ A^tot to estimate 
n(H2) for an assumed depth along the line of sight. At 
the adopted distance of NGC 2023, the SR angular width 
of ~2" corresponds to ~10"'^^ cm, which closely approxi- 
mates the model depth of the warm emitting region. Em- 
ploying /e « 5 we have ^(Ha) w A^to°t''/5 x 10^^ « 5 x 10"^ 
cm~^. The corresponding nu = 2 x ^(Ha) is, of course, 
higher. This estimate, which is PDR model indepen- 
dent, provides a strong support for the case that the SR 
is a relatively high-density region, and that our PDR 
modeling computes very reasonable values for the two 
parameters nn and x- Note that, just like Nq°^'^, the ob- 
servational estimates of both iV™/'^ and n(Il2) are lower 
limits owing to the underestimation of the unobserved 
(but higher-valued) J — and 1 level populations. 

3.4. Atomic Line Emission from Si and Fe 

In § 2.3 we remarked that the only two PDR-generated 
atomic emission lines present in our NGC 2023 spectra 
belong to [Si II] and [Fe II]. The former fine at 34.82 ^m 
is detected from the SR by Spitzer with an average inten- 
sity of 2.1 X 10~^ erg s~^ cm~^ sr~^. This is tw i ce the 
observed intensity estimated by IKaufman et al.l (|2OO60 
based on the ISO data,. 

The Kaufman et aljl (I2OO60 model prediction for [Si II] 
intensity was 6 x 10~^ erg cm~^ sr~^ based on nn 
and X from ISO H2 observations. Their model used a 
gas phase silicon abundance of Si/H ^ 1.7 x 10^^ for a 
depletion of ^20 relative to the solar abundance. The 
difference between model and observed intensity was at- 
tributed to additional depletion by a factor of ^6 com- 
pared to the model value. There is a strong dependence 
of model line intensities on both nu and x so that our 
PDR parameter values, which are based on the Spitzer 
H2 observations (§ 3.2.2) lead to e ven higher predicted in- 
tensity values. Figures 1 and 2 of IKaufman et al.l (|2006D 
present (nn, x) grids for normal models, predicting [Si II] 
and [Fe II] emission intensities, respectively. We use 
those grids and depletions for comparison with observed 
values, and do not repeat such calculations here since, as 
far as atomic emission computations ar e concerned, es - 
sentially the same models used by .Kaufman et al.l (|2006( ) 



are used here. 

Our H2 mapping of the SR indicated a log (nn, X.) solu- 
tion of (5.2, 4.0), for which the predicted [Si II] intensity 
is ^10"'^ erg s~^ cm"^ sr^^. Next, this value needs to be 
corrected for beam filling. As a rough estima te we use the 
ratio of the diffraction limit at 35 ^ui (8'.^5. iHouck et al] 
I2004f) over the SR limit of 2" to get -2 x 10"'* erg s'^ 
cm~^ sr~^. Consequently, [Si II] predictions based on 
our PDR model maps are '^10 times higher than the ob- 
served value. To bring our observations and model into 
agreement, the gas-phase Si/H abundance would have 
to be '^1.7 X 10"'', i.e., depleted by a factor of ^200 
relative to the solar Si abundance. Our results arc con- 
sistent with pre vious investigations (Draine & Bertold^ 
I2OOOI : iKaufmaiTet al.l l2006l ) in finding that a large de- 
pletion of gas-phase Si is required to bring models and 
observations into agreement. 

As for the weaker [Fe II] line at 25.99 /xm, it has 25% of 
the emission intensity of [Si II], or 5 x 10~^ erg s~^ cm~^ 
sr^^. This is the first reported detection of [Fe II] emis- 
sion in NGC 2023. The model-predicted [Fe II] to [Si II] 
intensity ratio is ^0.1, indicating that the Fe/Si abun- 
dance ratio is ~2.5 times higher than the values adopted 
in the models. Since the Si requires further depletion 
by '^10 we expect the Fe requires further depletion by 
~10/2.5 = 4. Indeed, the observed [Fe II] line inten- 
sity is well below the predicted value based on the model 
abundance of Fe/H = 1.7x10"^, and requires a gas phase 
abundance of Fe/H 4 x 10^^, or ^1/800 times the so- 
lar value. The depletion is ^4 — 6 times higher than that 
observed in diffuse clouds (jSavage &: Sembach 1996 ) and 
^,50 times higher t han that estimated in several PDRs 
(jOkada et al.ll2008l) . Our anomalously high Fe depletion 
could be a reflection of the higher value of on-SR gas 
density relative to densities probed in previous studies. 

4. GAS HEATING PROCESSES 

As was mentioned in the Introduction, previous studies 
have indicated the presence of low- J Tox values that seem 
to be rather high and thus pose a challenge to mod els of 
gas heating in PDR environments. .Tinimernian n et al.l 
U99g) found Tex - 500 K in S140 and were able to model 
observed line intensities by employing an initial gas tem- 
perature (at the computational edge of the cloud) of Tn 
= 1000 K, as well as cos(6') = 0.1. iFuente et all ([19991 ) 
found a range of T^^ « 300—700 K in NGC 7023. T oward 
the southern hah of NGC 2023, iFleming et all (poTol ) 
found Tex ~ 500 — 1400 by fitting H2 excitation curves 
with a single (hot) component ove r an area larg e r than 
the one considered here. Finally, 'Haba rt et al.l (|201lD 
analyzed Spitzer data from a few PDRs with low-valued 
FUV fields and found OPR-independent (AJ = 2) Tex 
values between 200 — 750 K toward the northern half of 
NGC 2023. 

Our data toward NGC 2023 provide de-zigzagged val- 
ues of Tox that range over 240 — 700 K for on-SR readings 
and essentially overlap all the pure-rotational results for 
other PDRs mentioned above. Furthermore, our PDR 
modeling, which consistently solves for the temperature 
structure of the gas for each PDR depth layer, shows a 
Tgas range of 300 — 750 K over the H2 line formation re- 
gion of, e.g., Av ^ 0.5 — 2. The close correspondence be- 
tween Tex and Tgas values calculated by the PDR models 
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confirms that tlie pure-rotational H2 emission lines de- 
tected and analyzed here are primarily thermally excited 
by collisions rather than radiatively excited by FUV pho- 
tons. 

Figure [T3] shows the depth variations of prominent gas 
heating processes for parameter values nn = 2 x 10^ cm^'^ 
and X = lO'*- These processes include grain photoelectric 
heating as well as heating via collisional de-excitation of 
FUV pumped H2, and via H2 dissociation. Clearly, grain 
photoelectric heating is dominant throughout the PDR 
layers where H2 line emission arises, with other heating 
processes contributing <25% of the heating budget. The 
successful reproduction of H2 data toward NGC 2023 im- 
plies that FUV interactions with the dust and gas com- 
ponents of this PDR do not require an additional energy 
input in the form of mechanical heating from, e.g., tur- 
bulence or shocks in this source. 

5. CONCLUDING REMARKS 

The rich molecular spectrum of H2 provides a rigor- 
ous diagnostic tool in the study of PDRs, the FUV- 
irradiated envelopes of molecular clouds. We showed 
that very good agreement can be obtained between mod- 
eled and observed absolute values of the line intensities 
of rotationally-excited H2 toward the SR in NGC 2023. 
According to our PDR models, the highest values of H2 
emission, which emanate from the narrow SR, require 
densities up to ~2 x 10^ cm~^ and radiation fields up 
to ~10^ times the local Galactic field. These values are 
well within the observationally acceptable range and are 
consistent with other PDR observables. The agreement 
between data and models is a direct result of improved 
sampling of the emission owing to the fine spatial res- 
olution of Spitzer instruments, as well as of recent im- 
provements in PDR modeling, including a more detailed 
treatment of H2, and an enhanced H2 formation rate 
for PDRs. Our re sults do not confirm the finding by 
iHabart et all (|201in of order-of- magnitude discrepancies 
between PDR model results and Spitzer H2 data, which 
could be explained by their observations of PDRs illumi- 
nated by mainly low-FUV fields. In contrast, for NGC 
2023 the FUV field is sufficiently high to dominate any 
non radiative heating that might be present. 

The fact that our model gives a good match to H2 
rotational line intensities, and to their associated run 
of increasing Tqx with J, is a good indication that the 
model includes an adequate treatment of heating and 
cooling processes. In particular, the dominant heating 
process via photoelectrons is sufficient to maintain the 
correct Tgas profile and II2 emission distribution within 
the PDR, without additional mechanical sources of heat- 
ing. Furthermore, an OPR resulting from H collisions 
and grain accretion provides a value that is ^3 and is 
matched by fitting both observed and modeled H2 col- 
umn densities. Thus our steady state computation pro- 
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duces the observed ratios among even- and odd- J states, 
without the need to artificially adjust the OPR. 

According to the PDR model, collisional excitation 
dominates the pure-rotational emission lines studied 
here. Thus the curvature of rising Tex with J is not the 
result of FUV pumping but of collisional thermalization 
of H2 levels in tandem with rising gas temperature within 
shallower layers of the PDR. In effect, the small fraction 
of radiative energy from HD 37903 that photoelectrically 
heats the gas is more important in controlling rotational 
level populations via collisions than FUV fluorescence. 

Our IRS maps show that the decrease in II2 emission 
intensity away from the SR is accompanied by reductions 
in riH, X, and fee, the latter being the ratio of data to 
face-on, beam-filled, normally-illuminated model inten- 
sity. The analysis of the four factors that affect the value 
of /off indicated that its reduction may be facilitated by 
a combined reduction in inclination and incidence fac- 
tors, which are countered by a more modest increase in 
the value of the beam-filling factor. Despite the difficulty 
of constraining all factors, the ob servations were s hown 
to be consistent with a previously (jField et al.lll994l) sug- 
gested 3-D structure of the region around the SR, namely 
a quasi-pyramidal molecular cloud towering above the 
FUV-carved cavity toward HD 37903. 

The successful interplay between PDR observations 
and theory supports the description of NGC 2023 as a par 
excellence example of a photodissociation region. This is 
especially significant because the results presented here 
require the shorter distance scale for HD 37903. Indeed, 
a very recent re-examination of the SIMBAD database 
on 2011 May 20 at 12:20 EDT revealed a newly pub- 
lished revision of th e HIPPARCOS parallaxes based on 
Ivan LeeuwenI (|2007[ ) . One of these parallaxes provides a 
new distance to HD 37903 of 300lg|!)'' pc, predicting an 

on-SR X of (l-03to'48) x lO'^, which is in extremely good 
agreement with our PDR modeling results. 
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TABLE 1 

Log of Spitzer Observations of NGC 2023 



Module 


Pixel 


Slit 


AOR 


Object 


Date 


a(J2000) 


<5(J2000) 


Exposures 




(") 


(" X ") 








(# X s) 


SH 


2.26 


4.7 X 11.3 


14033920 


NGC 2023 


2006/03/19 


05:41:37.63 


-02:16:42.6 


144 X 30 


LH 


4.46 


11.1 X 22.3 


14034176 


NGC 2023 


2006/03/19 


05:41:37.63 


-02:16:42.6 


30 X 60 


SL 


1.85 


3.7 X 57 


17977856 


NGC 2023 


2007/10/08 


05:41:37.63 


-02:16:42.6 


54 X 28 


SL 


1.85 


3.7 X 57 


17978112 


Sky bkgd 


2007/10/08 


05:40:26.21 


-02:54:40.3 


4 X 28 



TABLE 2 

GAUSSFIT Parameters'' and their Uncertainties*' 



Line 


Arest 




Ai 


Al — Arest 


A2 








{urn) 


(10-5) 


(ixm) 


(/im) 


(//m) 




(10-5) 


S(0) 








Map Median (off-SR) 






28.2188 


1.35(8) 


28.22.5(1) 


+0.006 


0.015(1) 


804(65) 


2.3(2) 


S{1) 


17.0348 


5.84(9) 


17.0378(2) 


+0.003 


0.0106(2) 


680(12) 


11.6(3) 


S{2) 


12.2786 


6.9(1) 


12.2801(1) 


+0.001 


0.0071(1) 


736(14) 


13.5(4) 


S{3) 


9.6649 


9.2(2) 


9.663(1) 


-0.002 


0.053(2) 


78(2) 


19.6(8) 


S(5) 


6.9095 


4.1(5) 


6.905(3) 


-0.005 


0.030(4) 


98(12) 


9(2) 


R.u." 




2—12% 


0.001—0.04% 




2—12% 


2—12% 


3—18% 



S(0) 
S(l) 
S(2) 
S(3) 
S(5) 



28.2188 
17.0348 
12.2786 
9.6649 
6.9095 



1-7(1) 
14.8(2) 
18.8(2) 
40.8(4) 
22(1)_ 
1—6% 



28.2231(7) 
17.0380(2) 
12.2798(1) 
9.6689(5) 
6.912(2) 
0.001—0.02% 



Single Pixel [7:8] (on-SR) 
+0.004 



+0.003 
+0.0008 
+0.004 
+0.002 



0.014(1) 
0.0105(2) 
0.0068(1) 
0.0511(6) 
0.031(2) 
1—8% 



844(64) 
688(13) 
772(13) 
80.3(9) 
95(6) 
1—8% 



2.8(3) 
29.4(7) 
33.8(7) 
84(1) 
54(4) 
1—9% 



R.U.' 



Each line is fitted with /(A) = Aaexp(—0.5[X — AiJ^/Aj) + A3 + A4A; the continuum (last two terms) is subtracted from the fit. 
Uncertainties for the last digits are in parentheses. 

Units for Aq are erg s'^ cm"^ sr-^; SAq = integrated line intensity, with conservative uncertainty taken from Aq & A2 in quadratures. 
R= Arest/(2.355A2) is the spectral resolution. 
" R.u. = Relative uncertainty. 



TABLE 3 

Normal Model Mapping of NGC 2023-South 



Map region 




X 


RMSD 


Jeff 


X/"H 


(# of pixels) 


(dex) 




(dcx) 




Global (170) 






Model: log nn = 5.0; log x = 3.7 






5.0 


5000. 


0.104 


0.50 


0.05 


On-SR (15) 


5.0 


5000. 


0.144 


1.07 


0.05 


Off-SR-S (9) 


5.0 


5000. 


0.104 


0.60 


0.05 


Off-SR-N (12) 


5.0 


5000. 


0.134 


0.37 


0.05 


Global (170) 






Grid: log nn = 4.7—5.3; log x = 3-6- 


-4.2 




5.0 


8000. 


0.062 


0.49 


0.08 


On-SR (15) 


5.2 


10000. 


0.062 


0.74 


0.06 


Off-SR-S (9) 


5.1 


6000. 


0.075 


0.48 


0.05 


Off-SR-N (12) 


4.9 


6000. 


0.064 


0.49 


0.08 



Note. — Values are sample medians. 
Number density units are cm-^. 
^ Ratio of data/model. 
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Fig. 1. — Left panel shows a near-infrared view of the southern half of NGC 2023 as combined from J ("blue"), H ("green"), and K 
("red") exposures. This highly magnified region of size 3' X 3' shows HD 37903 (Sellgren A) as the brightest star toward the top, and 
a pink-colored Southern Ridge (SR) just below the center of the image. Additional stars are identified by their Sellgren letters. Credit: 
ESO/J. Emerson/VISTA/Cambridgc Astronomical Survey Unit (ESO release 0949). Right panel shows a magnified region of 30" X 30" 
from an HST-ACS image, providing the highest angular resolution view of the SR to date. A blue square shows the size of an LH pixel 
from Spitzer-TRS. All figure labels were inserted manually and should not be presumed to have a level of precision better than 10%. Credit: 
|http: / / en. wikipedia.org/wiki/NGC_2023| based on ACS data set jSmwOl. 
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Fig. 2. — IRS all-module common field (irregular outline) overlaying the IRAC channel 4 image of NGC 2023, which is dominated by 
PAH emission at 8 /im. All maps to follow shall employ the instrumental orientation (vectors X, Y). As in left panel of Figure [T] Sellgren 
A is the IR-faint HE) 37903, whereas Sellgren C and D are two IR-bright young stellar objects. S is the SR, a narrow H2 emission filament, 
whereas SSE is the South-Southeastern Ridge, a wider emission clump. 
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Fig. 3. — Intensity maps of H2 line emission (uncorrected for extinction) toward the SR of NGC 2023. The color bars show the intensity 
scale in units of 10~^ erg s~^ cm"'^ sr~^. Each box spans 14 X 13 pixels, or ^ 1' X 1'. First five panels from top left show detected 
transitions from J = 2, 3, 4, 5 and 7 in the observed LH frame. Last panel employs the total intensity of all five emission lines as a 
background for contours of continuum intensity from the 8 fira IRAC image shown in Figure |2] The '+' indicates the position of Sellgren 
C. 
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Fig. 4. — Full spectral coverage from the four IRS modules SL2, SLl, SH, and LH toward NGC 2023, as obtained by averaging 15 pixels 
that sample SR emission. H2 and atomic emission lines are identified, and their rest wavelengths are indicated by dotted lines. Strong 
PAH features and dust continuum are evident. 
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Fig. 5. — Observed Nj values (symbols) for the single on-SR pixel LH[7:8] compared with model results (lines). Boxes show ANj = 
±20% corrected for Ay = 4.0 mag and triangles show the data prior to extinction correction. The dashed lino shows the unshifted log rajj = 
5.3 and x = 10* model. The solid line shows the same model shifted by —0.12 dex, or /^ff = 0.76 (see text), following R,MSD minimization. 
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Fig. 6. — Mapping of NGC 2023 with a grid of normal models. Upper left: logarithmic map of gas density variations. Upper right: 
linear map of FUV flux variations. Lower left: RMSD values are <0.16 dex over the map. Lower right: ratio of data to model required to 
minimize the RMSD. Each box spans 14 X 13 pixels, or 1' X 1'. Contours show the 30, 50 and 75% levels of the total intensity of all five 
H2 emission lines (from last panel of Figure[3ll. Red arrows extend exactly one half the distance toward HD 37903. 




Fig. 7. — Product fgf^ = cos(0)/cos(9) is contoured as a function of a and 6. As depicted in the lower right corner, a is the angle at 
the source between HD 37903 (*) and our line of sight (ffi), 6 is the surface inclination of the source, whose normal is indicated by and 
tp is the angle of FUV incidence. All values above upper limit of color bar (black contour of fgf,f, = 5) are mono-colored red. 




Fig. 8. — Cross-sectional cartoon of the suggested 3-D structure of NGC 2023 in the vicinity of the SR. The SR sits on top of a quasi- 
pyramidal high-density molecular core, the surface of which is made up of parallel cloud ridges that are subject to reduced levels of FUV 
irradiation. This core harbors the heavily obscured formation site of Sellgren D, see Figures [l] and |2] 
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Fig. 9, — Excitation diagram for the on-SR pixel LH[7:8]. Data boxes are ±20% in vertical extent, owing to IRS uncertainty. The data 
were de-zigzagged by shifting the odd-J values for OPR = 3 (triangles) toward the even-J levels, while minimizing the RMSD between the 
smooth excitation curve and the data. This pixel is found to have OPR ~ 1.9, with "cold" and "hot" Tex of 206 and 685 K, respectively 
(dashed lines), see also Figure [T2l 




Fig. 10.— Spatial variation of the OPR over the LH field toward NGC 2023. Box size is 14 X 13 pixels, or ?s 1' X 1'. Total H2 on-SR 
intensity in indicated by contour levels of 30, 50 and 75%, with the red arrow extending half the distance toward HD 37903 and the location 
of star C marked by '-|-'. 
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Fig. 11. — Three OPR-controUing processes and their sum are plotted as a function of depth into the PDR. Each process is found to 
dominate 0-H2 ^ P-H2 conversions over a certain range of Ay- At Ay < 1-5, P-H2 — >■ 0-H2 domination maintains an OPR of ft!3 (see 
scale on the right), but deeper into the PDR, 0-H2 — > P-H2 dominates owing to H2 accretion onto cold dust and is responsible for a rapid 
decline of the OPR. A horizontal dashed line shows the level of OPR = 1. 
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Fig. 13. — Modeled heating processes in the gas for a normal model with njj = 2 X 10^ cm~'^ and x = 10*- Heating is controlled by 
photoelectric emission from grains ("PE", solid curve) throughout the H2 line formation region indicated by colored emissivity (cooling) 
curves. Lesser heating contributions are from H2 dissociation (dotted curve) and from H2 ro- vibrational de-excitation (dashed curve plotted 
is net heating minus cooling). The dot-dashed curve shows the gas temperature profile, with scale on the right. 



